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We consider the rheology of soft-core frictionless disks in two dimensions in the neighborhood of 
the athermal jamming transition. From numerical simulations of bidisperse, overdamped, particles, 
we argue that the divergence of the viscosity below jamming is characteristic of the hard-core limit, 
independent of the particular soft-core interaction. We develop a mapping from soft-core to hard- 
core particles that recovers all the critical behavior found in earlier scaling analyses. Using this 
mapping we derive a relation that gives the exponent of the non-linear Herschel-Bulkley rheology 
above jamming in terms of the exponent of the diverging viscosity below jamming. 
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A variety of disordered soft solids, such as foams, col- 
loids and emulsions, are empirically observed to obey a 
non-linear rheology under steady state shearing known 
as the Herschel-Bulkley (HB) law pQ, a = cfq + cj b . Here 
7 is the shear strain rate, a is the average shear stress, cjq 
is the yield stress at 7 — >> 0, and b is an exponent usually 
(but not always) found experimentally to be in the range 
0.33 to 0.5 [2 . Detailed microscopic models of the viscous 
interaction in foams and emulsions have been studied to 
try and understand the HB form [3 . However in the 
limit of very slow strain rates, 7 —> 0, it seems likely that 
the rheology will be determined by collective effects and 
may be characterized as a critical phenomenon [3H6]. In 
this limit, the HB rheology has been treated in terms of 
a phenomenological model of slow glassy relaxation [4], 
and more recently in terms of the nucleation of localized 
plastic events [5]. 

Here we investigate this problem using numerical simu- 
lations of a model of athermally sheared, frictionless, soft- 
core disks. Such systems display a sharp jamming tran- 
sition as the packing fraction qb increases. For <j) < 0j, 
the system is liquid-like: at sufficiently small 7 the rhe- 
ology is linear [7 with a finite shear viscosity n = a/j 
that diverges as n ~ \<p — as <p — >> <pj. For <p > <pj 

rheology is non-linear with a finite yield stress cfq. By 
numerically establishing a mapping from sheared soft- 
core particles to sheared hard-core particles, we propose 
a relation between the exponent b of the non-linear HB 
rheology above <pj and the exponent f3 of the diverging 
viscosity of the linear rheology below 

Our model [8: is one of N bidisperse soft-core disks 
in two dimensions (2D), with equal numbers of particles 
with radii ratio 1.4. The soft-core interaction between 
two overlapping particles i and j is V(r^) = (e/a)^-, 
where 5{j = (1 — Vij/dij) is the relative particle overlap. 
Here is the particles center to center distance, is 
the sum of their radii, and a = 2 or 5/2 for harmonic or 
Hertzian interaction, respectively. Lengths are measured 
in units of the small particle diameter d si and energy in 



units of e. We use Durian's "mean-field" dynamics [9 
of overdamped particles with a viscous dissipation with 
respect to the imposed average linear shear velocity flow, 

3 

Time is measured in units of d s /Ce. Lees-Edwards 
boundary conditions [10 induce a uniform shear strain 
rate 7. We use N = 65536 particles so that finite size 
effects are negligible. Simulating at fixed 7 and (/>, we 
compute the steady state time average of the elastic part 
of the pressure tensor j8] , to define the scalar pressure p 
and shear stress a. We consider the pressure analog of 
the shear viscosity r] p = p/7, rather than 77, and restrict 
our analysis to a very narrow range about (/>j, specifi- 
cally 0.835 < (j) < 0.846 and 7 < 10 -6 for harmonic, and 
7 < 10 -7 for Hertzian, so as to allow us to ignore effects 
due to corrections to scaling [TTJ [12] . 

First we demonstrate the existence of the hard core 
limit below qbj. In Fig. [T] we compare r] p for both har- 
monic and Hertzian interations, for small 7 in the linear 
rheology region. We see excellent agreement, showing 
that the 7^0 limit of n p is independent of the partic- 
ular soft-core interaction. For the strict hard-core limit, 
one expects that particles at different strain rates 7 fol- 
low the same path through phase space, only with dif- 
ferent velocities, oc 7. For overdamped particles this 
implies that the contact forces also obey oc 7, and 
hence p oc 7. One may think of p/j in athermal shear 
driven flow as analogous to the virial p/T of equilibrium 
hard-core particles. In Fig. [TJd we replot rj p vs — 0, 
using = 0.8433 as previously determined [TTJ [13]. We 
see a clear algebraic divergence of n p over four decades as 
(j) — » </>j, demonstrating that the exponent /3 is character- 
istic of the hard-core limit, independent of the particular 
soft-core interaction. 

We next consider behavior outside the linear rheology 
(hard core) region, showing that one can map soft-core 
particles at a given (j) and 7 onto an equivalent hard-core 
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FIG. 1: (color online) (a) r\ v = p/j vs for data below 0j in 
the linear rheology region. Open symbols are for the harmonic 
interaction, solid symbols are for Hertzian, (b) Same data as 
in (a) but plotted vs <j)j — <f> with <j)j — 0.8433. 



system at a lower e ff(0,7), i- e - Vpifaj) = ^(</>eff,7 ^ 
0) = ^ d ((/5 e ff), using a simple form for e fj. If this map- 
ping holds, then even outside the linear rheology region 
f] p will have the scaling form, 



Vp(<t>,j) = Vp d (M = A((j)j - ^ eff )- 



(2) 



with ^ d ((/>) given by the curve in FigJTJ Such a map- 
ping was suggested by early work p~4] I15j in equilibrium. 
More recently Berthier and Witten [16] combined such a 
e ff approach with critical scaling to study the equilib- 
rium glass transition of soft spheres. A related analysis 
was done by Xu et al. [17], while more recent works 
have sought to extend this mapping over wider ranges of 
pressure [18] and to systems with applied uniform shear 
strain rate 7 [19], though still at finite T. 

Here we apply these ideas to an athermal shear-driven 
system. We follow Berthier and Witten [16] and assume 
that e ff is set by the extent of particle overlaps, as mea- 
sured by the average interaction energy per particle E. 
We thus make the Ansatz, 



(3) 



with h(E) an appropriate function to be determined. 
Since E = when there are no overlaps, h(0) = 0. 

We can now determine h(E) asymptotically close to 
(j)j by applying two simple conditions on <j) e R : 



, 7 0) = , for (j) < 



(4) 



Since E — >• as 7 — >• below </5j, overlaps vanish and 
6. The second condition is, 



^eff 



for (/> > 



(5) 



Since p — > po > as 7 — > for all > 0j, then r] p ^ 00 
everywhere along the dynamic yield stress curve. In a 
hard-core system, r] p — » 00 only at = </>j (0 > <fij being 
excluded by the non-overlapping condition). Thus, as 



7 — > 0, all > 0j in a soft-core system must map onto 
^eff = <j>j of the equivalent hard-core system. 

If we similarly define Eo(<p) = £"(0, 7 — > 0), then 
Eqs. ([3| and Q imply h(E (</))) = <j> - <j>j for > 0j. 
Close to 0j, £0 scales to zero algebraically, Eo ~ (0 — 
We thus conclude that h(E) ~ E 1 ^ yE , and so, 



eff (0,7) = 0-c[^,7)] 1/2/£ 



(6) 



We test this mapping by measuring r] p and E at var- 
ious (j) and 7, and fitting our data to Eqs. ([2| and 
taking A, </5j, /3, c and de as free fitting parameters. In 
Fig. [2^i we show our raw data rj p vs (j) for the harmonic 
interaction, including points both above and below <pj] 
f] p decreases with increasing 7, showing that our data in- 
clude points outside the linear rheology region. In Fig.[2]3 
we show the results of our fit to the (/> e ff model, plot- 
ting f] p vs <j)j — e ff • We find an excellent data collapse, 
yielding 0j = 0.84328 db 0.00007, (3 = 2.58 ± 0.10 and 
Ue = 2.18 ±0.02 [12]. We can compare this to our earlier 
results [11] from a more general critical scaling analysis 
that was independent of any <j) e R assumption. Defining 
the pressure exponent y p by po(4>) ~ {(j) — <j)j) Vp , our ear- 
lier results gave, <j>j = 0.84347 ± 0.00020, /3 = vz - y p = 
2.77 ± 0.20 [20]; taking he = 2y P for harmonic interac- 
tion, we get He = 2.16 ± 0.06. The excellent agreement 
between the two approaches establishes the validity of 
our soft to hard-core mapping, <j) e R , for the range of data 
we simulate. 
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FIG. 2: (color online) (a) ?7 P = p/j vs for data above and 
below 0j for harmonic soft core particles, (b) Same data as 
in (a) but plotted vs <j)j — e ff- Dashed line in (a) and solid 
line in (b) is the fit to the model of Eqs. |2]) and (pi. Symbols 
in (b) correspond to the legend given in (a). 



Fig. [3] shows a similar analysis for the Hertzian inter- 
action. Outside the linear rheology region, the Hertzian 
T] p is smaller than for the harmonic case due to the 
softer repulsion of the Hertzian cores. Consequently, 
our Hertzian data generally lies further from the asymp- 
totic 7—^0 hard-core limit and thus is poorer at de- 
termining the critical behavior. However, since the pa- 
rameters A, <j)j and /3 defining 7?£ d in Eq. (||) are char- 
acteristic of the hard-core limit, only the parameters 
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c and yE defining <p e ^ in Eq. (|6| should vary as the 
soft-core interaction is changed. We therefore fix A, 
and P to the values found from the harmonic data, 
and allow only c and yE to vary. The fit, shown in 
Fig. [3}}, is excellent and gives yE — 2.70 ± 0.04. We find 
^Hertzian ^harmonic = 1.24 ±0.05, in good agreement with 
the ratio a Hertzian/ a harmonic = x 25. Since E is re lated 
to the average particle overlap 5 by E ~ S a , this obser- 
vation suggests S ~ (0 — 0j) 108 , common to all soft-core 
interactions. 
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FIG. 3: (color online) (a) r\ v = p/j vs <j> for data above and 
below cj)j for Hertzian soft core particles, (b) Same data as 
in (a) but plotted vs — e ff. Dashed line in (a) and solid 
line in (b) is the fit to the model of Eqs. |2| and |6j), using 
the same values of A, <j)j and /3 as found for the harmonic 
interaction model. Symbols in (b) correspond to the legend 
given in (a). 



We now discuss the implications of our <p e ^ mapping for 
the Herschel-Bulkley rheology above <j)j. For observables 
X such as p, a or E, the HB form for small 7 is, 



X^,j) = X o (0) + C x j b 



(7) 



where bx is the HB exponent, and Xq((/)) = X(4>, 7 — >• 
0). First we review some conclusions that follow from a 
general critical scaling ansatz [IT] , independent of our 
(j) e fi mapping. X is expected to have the scaling form, 



x(0,7) = i^r^± 



- 0J , (8) 



where zv = f3 + y p [20] and X± are the scaling functions 
for Scj) ^ 0. As — » 0j, and the argument of the scaling 
function diverges, the dependence of X on £0 should drop 
out, thus requiring, 

X±{x^kx>)~x vx I zv , (9) 

so that exactly at we have the non-linear rheology, 



X(0j,7) - j qx •> witn 4x 



zv 



yx 



. (10) 



Unlike /3, we see that the rheology exponent at 
does depend on the particular soft-core interaction, via 
the exponents y p and yx- 



For <j) > cj)j, as 7 — >> 0, Eq. ([7]) requires the scaling 
function to have the form, 

X + (x ^ 0) = c x + c' x x bx , (11) 

with cx-,c' x constants, so that we recover Eq. Q with 

X = c x H yx , and C x = c' x H yx ~ hxZV • (12) 

Thus the coefficient Cx of the HB law of Eq. ^ must 
have a particular scaling dependence on </> as <\> — » 0j. 

We now return to our </> e ff model and consider the pres- 
sure. From the definition of r\ v and Eq. Q we can write, 



p(0,7) 



7A 



7A 



(13) 



Substituting in Eq. ([7| for i£, and expanding h(E) to first 
order for small 7, we get, 



jA 



(14) 



t) to cancel 



where we used Ji(Eq) = cE { 
out the leading term in the expansion of h(E). As 7 — >• 
above 0j, p — » po is finite. We thus conclude from 
Eq. (14) that we must have, 6# = 



To determine the HB exponent for pressure, we just 



extend the expansion in Eq. (14) to next order, 



jA 



[h f (E )C E j bE + %h»(E )C*<y 2b B] 



Po 



_ f3h"(E )C E bE 
2h'(E ) 7 



(15) 
(16) 



Comparing to Eq. ^ for p we conclude that b p = bE = 
l/f3. Similar results hold for the yield stress a. We 
thus conclude that the HB exponents are all equal to 
b = 1//3 « 0.38 ± 0.02, and by our earlier arguments, 
they are all independent of the particular soft-core in- 
teraction. We note that a similar value, b ~ 0.36, was 
recently reported in experiments on sheared foams [2Tj . 

We next numerically check our prediction for b. In 
Fig. [4^i we show the scaling collapse of energy E according 
to Eq. ^ for the harmonic case. Using the parameters 
found from our e fT fit, we see an excellent data collapse. 
In Fig. 03 we plot 6 + = E /\5(j)\ yE -c E vs x = ^il\^\ zv for 
> 0j, where ce of Eq. ( pT| ) is obtained from c of Eq. (|6| 



via ce = l/c VE . From Eq. (|11| we expect £+ ~ x b at 
small x, while from Eq. (|9| we expect £+ ~ x gs at large 
x. We consider £^ rather than p since there is a greater 
difference between the exponents b and g# than between b 
and </p. Fitting the data of Fig.[4]3 separately at small and 
large x we find power-law behaviors with b w 0.37 and 
^ w 0.55, respectively, in reasonable agreement with the 
values expected from our e g analysis, b = 1//3 = 0.38 
and g# = yE/(/3 + y p ) = 0.59. The horizontal dashed line 
in Fig. [4b locates the value ce on the vertical axis. Data 
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FIG. 4: (color online) (a) Scaling collapse of energy as in 
Eq. ([§, using values 0j, and zv = /3 -\- y p obtained from 
our fit to the e ff model. Two branches of the scaling function 
correspond to <f> above and below <pj. (b) Plot of E/\5(j)\ yE — ce 
vs x = 7/|^(/)| 2zy for > 0j. Dashed line is a fit to the large x 
data, giving a power- law behavior with exponent qs ~ 0.55; 
solid line is a fit to the small x data, giving a power-law 
behavior with exponent b ~ 0.37. Horizontal short dashed line 
is the value ce] data below this line satisfy (E — Eq) / Eq < 1. 



below this line satisfy the condition (E — Eq) / Eq < 1. We 
see that this condition locates the crossover from small 
to large x behavior, which occurs near x ~ 10 3 , or when 
7 « 10 3 <^. 

Although we go to smaller values of 7 than are typi- 
cally used in other works, the closest our data for <p> <pj 
approaches the yield stress line is (p—po)/po ^ 0.18. One 
can always question whether this is close enough to give 
the true asymptotic critical behavior, or whether rheolog- 
ical behavior might change at even smaller 7. We leave 
further investigation of this point to future work. Here 
we note that experimental fits to the HB form usually 
involve data extending well above this limit down to val- 
ues that typical do not go below ~ 0.1 [22] . Thus, even if 
our (/> e ff model ultimately breaks down closer to the yield 
stress line, our results remain of considerable relevance 
for understanding the experimentally determined value 
of the HB exponent in numerous physical systems. 

Our analysis has been for a model with dissipation to 
an external reservoir, yielding a Newtonian (linear) rhe- 
ology below <j)j. In athermal granular systems, with col- 
lisional dissipation and inertial effects, one expects Bag- 
nold scaling [23], <r,p ~ j 2 . In this case the Bagnold 
coefficient scales as B p = p/j 2 ~ (<pj — </>) _/3 as jam- 
ming is approached [24]. If a similar <j) e ^ model holds, 
one can repeat all the steps of our above argument to 
arrive at the HB exponent for this case, b = 2/ ' j3' ', while 
exactly at <pj we have qx = 2yx / {Vp+P')- From Ref. [24] 
we expect f3' w 4, giving b w 0.5. We note that the expo- 
nent b w 0.5 was observed in recent molecular dynamic 
simulations of a 2D athermal Lenard- Jones glass [5]. 

To conclude, by mapping soft-core particles at general 
(0,7) to hard core particles at (</> e ff,7 — > 0), we map the 
nonlinear rheology as 7 — >> above jamming to the linear 



rheology as 7 —> below jamming, resulting in a relation 
between the HB exponent b and the viscosity exponent 
/3. When comparing our results to experiments, how- 
ever, several points must be kept in mind: (i) The HB 
exponent b found here characterizes the rheology only for 
sufficiently small 7. Near 0j, as 7 increases, one crosses 



into a region characterized by the exponent q of Eq. (10) 
(see Fig. [4|. For systems with significant collisional dissi- 
pation and inertial effects, a crossover from Newtonian to 
Bagnold rheology is also possible [25H27] . Fitting the HB 
form to data that spans such crossover regions will there- 
fore result in an effective exponent b different from that 
reported here, (ii) The numerical value of b we report 
here results from the simple Durian "mean-field" model 
of dissipation, Eq. ([!]). Different models for viscous dissi- 
pation may yield different values for the exponent /?, and 
hence for b [23 [28]. 
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